NeXL is a collection of packages (libraries) for the Julia language that implement common X-ray microanalysis algorithms, physical data and utilities.
The core NeXL packages are
NeXL builds on many packages from the Julia infrastructure including
Julia is a high-performance scripting language designed for numerical data analysis. The base Julia libraries provide basic scalar, vector and matrix numerical analysis tools. Julia differs from most scriptable languages in that all code is compiled before execution. This leads to a "slow then fast" user experience. The first time a method is used, it and all the methods it depends upon must be compiled. This can lead to a significant delay. However, once compiled, the code is fast - often as fast as C or Fortran code. You will notice this "slow then fast" nature as you work with this notebook. It is particularly evident when creating plots.
We will start by exploring NeXLCore since this package provides core functionality for the other NeXL packages.
using DrWatson
@quickactivate "IUMAS"
using NeXLCore
Activating project at `c:\Users\nritchie\Desktop\IUMAS`
We will first construct a datum representing an element and then manipulate this object to extract elemental data.
Julia is not an object-oriented language. Instead it uses what is called "multiple dispatch" to determine which method implementation to call dependent on the arguments handed to the function.
el = n"Tc" # Here we are using a special syntactic sugar to convert the string "Tc" into an datum
| category | transition metal |
|---|---|
| atomic mass | 98.0 u |
| density | 11.0 g/cm³ |
| molar heat | 24.27 J/mol⋅K |
| melting point | 2430.0 K |
| boiling point | 4538.0 K |
| phase | Solid |
| shells | [2, 8, 18, 13, 2] |
| electron configuration | 1s² 2s² 2p⁶ 3s² 3p⁶ 4s² 3d¹⁰ 4p⁶ 5s² 4d⁵ |
| appearance | shiny gray metal |
| summary | Technetium (/tɛkˈniːʃiəm/) is a chemical element with symbol Tc and atomic number 43. It is the element with the lowest atomic number in the periodic table that has no stable isotopes:every form of it is radioactive. Nearly all technetium is produced synthetically, and only minute amounts are found in nature. |
| discovered by | Emilio Segrè |
| source | https://en.wikipedia.org/wiki/Technetium |
Let's get some properties of this elemental datum.
z(el), a(el), symbol(el), name(el)
(43, 98.0, "Tc", "Technetium")
We can combine elements into materials. The mat"..." macro converts expressions into Material structs.
mat = mat"Ca(PO4)3OH"
mat, typeof(mat)
(Ca(PO4)3OH[P=0.2717,Ca=0.1172,O=0.6082,H=0.0029], Material{Float64, Float64})
An alternative syntax allows you to specify additional properties of the Material.
mat = parse(Material, "Ca(PO4)3OH", name = "Apatite", density = 2.2)
Apatite[P=0.2717,Ca=0.1172,O=0.6082,H=0.0029,2.20 g/cm³]
In addition to being able to parse chemical formulae, NeXLCore can also parse expressions like the next one. In this case, the multiplier represents the mass fraction of the element. The right hand side can also be a chemical formula as though the material was made from various mass fractions of consituent materials.
adm = parse(Material, "0.3399*O+0.0664*Al+0.0405*Si+0.0683*Ca+0.0713*Ti+0.1055*Zn+0.3037*Ge", name="ADM-6005a")
ADM-6005a[Ti=0.0713,Al=0.0664,Ge=0.3037,Ca=0.0683,Zn=0.1055,O=0.3399,Si=0.0405]
In addition to elements and materials, NeXLCore provides mechanisms to work with characteristic X-rays, atomic shells and other atomic physics data. The n"..." macro can be used to create Element, CharXRay, AtomicSubShell and SubShell structs.
cxr = n"Tc K-L3"
cxr, typeof(cxr)
(Tc K-L3, CharXRay)
ass = n"Tc K"
ass, typeof(ass)
(Tc K, AtomicSubShell)
Let's get some properties of the CharXRay datum. You'll notice that the energy is in electron-volts as is the case for all energies within NeXL.
energy(cxr), inner(cxr), energy(inner(cxr)), outer(cxr), energy(outer(cxr)), typeof(inner(cxr))
(18367.0, Tc K, 21050.0, Tc L3, 2683.0, AtomicSubShell)
To determine every characteristic X-ray generated by a particular element, you can use the characteristic(...). To generate sub-sets of the transitions replace alltransitions with ktransitions, ltransitions or mtransitions.
cxrs=characteristic(el, alltransitions)
cxrs, length(cxrs)
(Tc K-L3 + 9 others, Tc L3-M5 + 23 others & Tc M5-N3 + 9 others, 44)
Now we can use various functions to extract X-ray related physical data from elements and materials. The mac(...) function returns the mass absorption coeffficient in $cm^2/g$.
mac(mat, n"Ca K-L3"), mac(n"Ca", n"Ca K-L3"), mac(mat, 3692.0), mac(n"Ca", 3692.0)
(257.51384578832244, 139.20100028283625, 257.51384578832244, 139.20100028283625)
This is an example of one of Julia's core design features - multiple dispatch. The function mac(...) is called on four distinct argument types.
mac(Material, CharXRay)mac(Element, CharXRay)mac(Material, Float64)mac(Element, Float64)This is similar but extends on the way object oriented code can use the object type to determine which function implementation to call.
NeXL uses the Gadfly library to plot spectra and other data items.
using Gadfly
plot(e->mac(mat, e), 100.0, 1.0e4, Scale.y_log10)
plot(
layer(e->mac(mat, e), 100.0, 1.0e4, Theme(default_color="Red")),
layer(e->mac(n"Ca", e), 100.0, 1.0e4, Theme(default_color="Green")),
layer(e->mac(n"P", e), 100.0, 1.0e4, Theme(default_color="Blue")),
Scale.y_log10
)
There is much more functionality in NeXLCore, but that is enought to get us started. Next let's take a look at NeXLSpectrum.
NeXLSpectrum defines methods for reading, writing, manipulating, plotting and fitting X-ray spectra and X-ray spectrum images ("hyperspectra").
using NeXLSpectrum
[ Info: Loading Gadfly support into NeXLSpectrum.
First, we will read a spectrum from disk using the loadspectrum(...) method. It take a filename and is able to sniff the file type for many common spectrum file formats.
sp = loadspectrum(joinpath(datadir(), "exp_raw", "ADM6005a spectra","ADM-6005a_1.msa"))
* 128860.0
*
*
*
*
*
* *
* *
* *
* ** **
* ** ***** ** **
**************************************************************************************************** 20.0 keV
Spectrum{Float64}[ADM-6005a_1, -484.20818 + 5.01716⋅ch eV, 4096 ch, 20.0 keV, unknown composition, 6.81e6 counts]
set_default_plot_size(8inch, 4inch)
specs=loadspectrum.(joinpath(datadir(), "exp_raw", "ADM6005a spectra","ADM-6005a_$i.msa") for i in 1:15)
plot(specs..., klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"])
plot(specs..., klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"], xmax=6.0e3)
To quantify this data, we need standard spectra. The standard spectra are preprocessed to optimize the fitting process. The spectrum files are read from disk and must contain the following properties :BeamEnergy, :LiveTime, :ProbeCurrent and :TakeOffAngle to be fit and quantified. It is possible to pre-load the standard spectra and edit the properties in the script if they are not present in the spectrum file.
path=joinpath(datadir(), "exp_raw", "ADM6005a spectra")
refs=references( [
reference(n"O", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
reference(n"Si", joinpath(path, "SiO2 std.msa"), mat"SiO2"),
reference(n"Al", joinpath(path, "Al std.msa"), mat"Al"),
reference(n"Ca", joinpath(path, "CaF2 std.msa"), mat"CaF2"),
reference(n"Ti", joinpath(path, "Ti trimmed.msa"), mat"Ti"),
reference(n"Zn", joinpath(path, "Zn std.msa"), mat"Zn"),
reference(n"Ge", joinpath(path, "Ge std.msa"), mat"Ge"),
], 132.0
)
using DataFrames
asa(DataFrame, refs)
| Row | Spectrum | Beam Energy (keV) | Probe Current (nA) | Live Time (s) | Material | Lines | ROI | P-to-B | S-to-N |
|---|---|---|---|---|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | Material… | Array… | UnitRang… | Float64 | Float64 | |
| 1 | SiO2 | 20.0 | 5.0 | 100.0 | SiO2[Si=0.4674,O=0.5326] | O K-L3 + 1 other | 175:228 | 21.8788 | 4119.94 |
| 2 | SiO2 | 20.0 | 5.0 | 100.0 | SiO2[Si=0.4674,O=0.5326] | Si K-L3 + 3 others | 409:491 | 20.0938 | 8705.35 |
| 3 | Al | 20.0 | 5.0 | 100.001 | Al[Al=1.0000] | Al K-L3 + 3 others | 360:433 | 19.8537 | 14014.8 |
| 4 | CaF2 | 20.0 | 2.5 | 100.0 | CaF2[Ca=0.5133,F=0.4867] | Ca K-L3 + 3 others | 788:938 | 19.4466 | 5767.9 |
| 5 | Ti | 20.0 | 5.0 | 100.001 | Ti[Ti=1.0000] | Ti L3-M5 + 13 others | 153:226 | 4.7057 | 1092.39 |
| 6 | Ti | 20.0 | 5.0 | 100.001 | Ti[Ti=1.0000] | Ti K-L3 + 3 others | 948:1125 | 22.533 | 10719.4 |
| 7 | Zn | 20.0 | 5.0 | 100.002 | Zn[Zn=1.0000] | Zn L3-M5 + 13 others | 247:348 | 17.6177 | 9149.2 |
| 8 | Zn | 20.0 | 5.0 | 100.002 | Zn[Zn=1.0000] | Zn K-L3 + 1 other | 1753:1882 | 13.61 | 4272.31 |
| 9 | Zn | 20.0 | 5.0 | 100.002 | Zn[Zn=1.0000] | Zn K-M3 + 3 others | 1945:2065 | 2.6177 | 686.219 |
| 10 | Ge | 20.0 | 5.0 | 100.001 | Ge[Ge=1.0000] | Ge L3-M5 + 15 others | 277:397 | 14.9005 | 9088.97 |
| 11 | Ge | 20.0 | 5.0 | 100.001 | Ge[Ge=1.0000] | Ge K-L3 + 1 other | 1997:2135 | 10.5729 | 3025.67 |
| 12 | Ge | 20.0 | 5.0 | 100.001 | Ge[Ge=1.0000] | Ge K-M3 + 5 others | 2223:2358 | 1.93878 | 478.809 |
The references are then used to fit the unknown spectra and matrix correction is performed. The quantify(...) method can be used to perform both operations or the fit_spectrum(...) method can be used to fit the references to the unknowns to produce k-ratios and then the k-ratios fed to the quantify(...) method to perform matrix correction. The later two-step process is useful when you want to access the k-ratio data.
qrs=map(sp->quantify(sp, refs),specs)
qdf=asa(DataFrame, qrs, nominal=adm)
| Row | Material | O | Al | Si | Ca | Ti | Zn | Ge | Total |
|---|---|---|---|---|---|---|---|---|---|
| String | Abstract… | Abstract… | Abstract… | Abstract… | Abstract… | Abstract… | Abstract… | Abstract… | |
| 1 | ADM-6005a_1 | 0.380759 | 0.0680223 | 0.0388277 | 0.0639499 | 0.0733595 | 0.114264 | 0.318852 | 1.05803 |
| 2 | ADM-6005a_2 | 0.379643 | 0.0678927 | 0.0389075 | 0.0639314 | 0.0728992 | 0.114022 | 0.319517 | 1.05681 |
| 3 | ADM-6005a_3 | 0.381752 | 0.0680904 | 0.0387023 | 0.0640476 | 0.0730505 | 0.114183 | 0.317511 | 1.05734 |
| 4 | ADM-6005a_4 | 0.38067 | 0.068038 | 0.0389573 | 0.0637997 | 0.073006 | 0.114154 | 0.31897 | 1.05759 |
| 5 | ADM-6005a_5 | 0.380851 | 0.0683536 | 0.0386646 | 0.0640086 | 0.0730713 | 0.114185 | 0.317727 | 1.05686 |
| 6 | ADM-6005a_6 | 0.381398 | 0.0682851 | 0.0388131 | 0.0636924 | 0.0729745 | 0.113707 | 0.317195 | 1.05607 |
| 7 | ADM-6005a_7 | 0.382109 | 0.0678588 | 0.0389113 | 0.0640267 | 0.0732265 | 0.113666 | 0.319045 | 1.05884 |
| 8 | ADM-6005a_8 | 0.381185 | 0.068217 | 0.0390383 | 0.064036 | 0.0731137 | 0.113816 | 0.319457 | 1.05886 |
| 9 | ADM-6005a_9 | 0.381249 | 0.0681775 | 0.0389637 | 0.0639798 | 0.0730421 | 0.113619 | 0.316589 | 1.05562 |
| 10 | ADM-6005a_10 | 0.380373 | 0.0683098 | 0.0389339 | 0.063795 | 0.0728326 | 0.114285 | 0.318074 | 1.0566 |
| 11 | ADM-6005a_11 | 0.381009 | 0.0680763 | 0.0390429 | 0.0639864 | 0.0729592 | 0.113997 | 0.318805 | 1.05788 |
| 12 | ADM-6005a_12 | 0.381774 | 0.067932 | 0.038931 | 0.0639285 | 0.0727576 | 0.113903 | 0.317832 | 1.05706 |
| 13 | ADM-6005a_13 | 0.381688 | 0.0681477 | 0.0388243 | 0.0640871 | 0.0729426 | 0.113803 | 0.318703 | 1.0582 |
| 14 | ADM-6005a_14 | 0.38168 | 0.0678288 | 0.0389896 | 0.0639993 | 0.0728709 | 0.113778 | 0.318276 | 1.05742 |
| 15 | ADM-6005a_15 | 0.382318 | 0.0682955 | 0.0387906 | 0.0640459 | 0.0731022 | 0.114217 | 0.318665 | 1.05943 |
| 16 | ADM-6005a | 0.3399 | 0.0664 | 0.0405 | 0.0683 | 0.0713 | 0.1055 | 0.3037 | 0.9956 |
We can use the describe(...) method from DataFrames to generate summary statistics.
describe(qdf[1:end-1,2:end], :mean, :std, :min, :max)
| Row | variable | mean | std | min | max |
|---|---|---|---|---|---|
| Symbol | Float64 | Float64 | Float64 | Float64 | |
| 1 | O | 0.381231 | 0.000704508 | 0.379643 | 0.382318 |
| 2 | Al | 0.0681017 | 0.000172422 | 0.0678288 | 0.0683536 |
| 3 | Si | 0.0388866 | 0.000113223 | 0.0386646 | 0.0390429 |
| 4 | Ca | 0.0639543 | 0.000110957 | 0.0636924 | 0.0640871 |
| 5 | Ti | 0.0730139 | 0.000153452 | 0.0727576 | 0.0733595 |
| 6 | Zn | 0.113973 | 0.0002319 | 0.113619 | 0.114285 |
| 7 | Ge | 0.318348 | 0.000845766 | 0.316589 | 0.319517 |
| 8 | Total | 1.05751 | 0.00106329 | 1.05562 | 1.05943 |
Or we can take a closer look at the data from one spectrum.
asa(DataFrame, qrs[1])
| Row | Label | Element | Standard | Xrays | C | ΔC | k | Δk | g | Z | A | F | c | gZAFc |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| String | String | String | String | Float64 | Float64 | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | Float64? | |
| 1 | ADM-6005a_1 | Ti | Ti | Ti K-L3 + 3 others | 0.07336 | 0.00012 | 0.0643276 | 0.000103282 | 1.0 | 0.925634 | 0.947328 | 1.0 | 1.0 | 0.876879 |
| 2 | ADM-6005a_1 | Al | Al | Al K-L3 + 3 others | 0.068022 | 0.00012 | 0.0280766 | 4.99883e-5 | 1.0 | 1.03396 | 0.398685 | 1.0013 | 1.0 | 0.412758 |
| 3 | ADM-6005a_1 | Ge | Ge | Ge K-L3 + 1 other | 0.31885 | 0.00064 | 0.263479 | 0.000525419 | 1.0 | 0.831431 | 0.993868 | 1.0 | 1.0 | 0.826333 |
| 4 | ADM-6005a_1 | Ca | CaF2 | Ca K-L3 + 3 others | 0.06395 | 8.9e-5 | 0.121389 | 0.000168661 | 1.0 | 1.04312 | 0.926501 | 1.00822 | 1.0 | 0.974395 |
| 5 | ADM-6005a_1 | Zn | Zn | Zn K-L3 + 1 other | 0.11426 | 0.00028 | 0.111712 | 0.000272437 | 1.0 | 0.879367 | 0.999501 | 1.11233 | 1.0 | 0.977661 |
| 6 | ADM-6005a_1 | Si | SiO2 | Si K-L3 + 3 others | 0.038828 | 9.1e-5 | 0.0530552 | 0.000124176 | 1.0 | 1.11012 | 0.575032 | 1.00057 | 1.0 | 0.63872 |
| 7 | ADM-6005a_1 | O | SiO2 | O K-L3 + 1 other | 0.38076 | 0.00053 | 0.487558 | 0.000672807 | 1.0 | 1.10908 | 0.614952 | 0.999868 | 1.0 | 0.681942 |
Associated with each spectrum are a collection of properties. These properties are often read from the spectrum file but can also be assigned manually.
NeXLCore.properties(sp)
Dict{Symbol, Any} with 11 entries:
:Owner => "Bruker Nano"
:RealTime => 112.108
:BeamEnergy => 20000.0
:Title => "Bruker Nano spectrum Acquisition"
:Elevation => 0.698132
:LiveTime => 100.001
:Filename => "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_raw\\ADM6…
:ProbeCurrent => 5.0
:TakeOffAngle => 0.698132
:Name => "ADM-6005a_1"
:AcquisitionTime => DateTime("2019-07-12T09:33:00")
Properties are identified with a Symbol which is represented by a colon ':' plus a name. There are many predefined spectrum properties but you can also define your own to associate custom data with spectra.
sp[:BeamEnergy], sp[:ProbeCurrent], sp[:LiveTime], sp[:TakeOffAngle], dose(sp)
(20000.0, 5.0, 100.001, 0.6981317007977318, 500.005)
This is one way to index spectra to extract data from the spectrum. It is also possible to extract channel count data using various different types as indices. To demonstrate this, we will use Integer, AbstractFloat or CharXRay types to index the same channel in the spectrum - the channel associated with the Ca K-L3 peak.
channel(n"Ca K-L3", sp), channel(Float64, n"Ca K-L3", sp), channel(energy(n"Ca K-L3"), sp), channel(Float64, energy(n"Ca K-L3"), sp)
(833, 833.384891053903, 833, 833.384891053903)
In the first case, we index the channel directly by channel index. In the second case, we use a floating point number which is assumed to represent an energy to index the same channel. Finally, we use a CharXRay to index the same channel.
sp[833], sp[energy(n"Ca K-L3")], sp[n"Ca K-L3"]
(18042.0, 18042.0, 18042.0)
Often it is useful to be able to determine which of a collection of spectra collected from a single material under similar conditions are similar to one another. The notion being that the 'best' spectra are those that are most similar. These can then be summed together to produce a single spectrum that best represents the material. The findsimilar(...) method extracts the similar spectra from a list and removes outliers.
plot(sum(findsimilar(specs, refs.detector, n"O")), klms=[n"O", n"Al", n"Si", n"Ca", n"Ti", n"Zn", n"Ge"])
The similarity is determined by comparing each spectrum with the sum of the other spectra and calculating a score that respects count statistics. Similar spectra have a similarity(...) score of approximately unity.
DataFrame(Name=name.(specs), S1=NeXLSpectrum.similarity(specs, refs.detector, adm), S2=NeXLSpectrum.similarity(specs), Counts=integrate.(specs))
| Row | Name | S1 | S2 | Counts |
|---|---|---|---|---|
| String | Float64 | Float64 | Float64 | |
| 1 | ADM-6005a_1 | 0.861481 | 1.04111 | 6.81189e6 |
| 2 | ADM-6005a_2 | 0.852062 | 0.970497 | 6.81126e6 |
| 3 | ADM-6005a_3 | 0.854675 | 0.99976 | 6.81006e6 |
| 4 | ADM-6005a_4 | 0.854278 | 0.980264 | 6.81625e6 |
| 5 | ADM-6005a_5 | 0.89041 | 0.980602 | 6.80921e6 |
| 6 | ADM-6005a_6 | 0.890928 | 0.962467 | 6.8086e6 |
| 7 | ADM-6005a_7 | 0.841213 | 0.976628 | 6.81929e6 |
| 8 | ADM-6005a_8 | 0.778826 | 0.971166 | 6.81616e6 |
| 9 | ADM-6005a_9 | 0.83867 | 1.01652 | 6.816e6 |
| 10 | ADM-6005a_10 | 0.847391 | 1.00353 | 6.81184e6 |
| 11 | ADM-6005a_11 | 0.861667 | 0.975445 | 6.81146e6 |
| 12 | ADM-6005a_12 | 0.823105 | 1.01421 | 6.81786e6 |
| 13 | ADM-6005a_13 | 0.814398 | 0.974283 | 6.81113e6 |
| 14 | ADM-6005a_14 | 0.777566 | 0.972348 | 6.81502e6 |
| 15 | ADM-6005a_15 | 0.780707 | 0.942442 | 6.81972e6 |
A HyperSpectrum represents a multi-dimensional array of point spectra. NeXLSpectrum can handle a arbitrary number of dimensions like 1-D (linescans), 2-D (area scans) or higher dimensions.
There are currently three ways to load a hyperspectrum.
We will use the DataDeps library to download the hyperspectrum data on demand from the NIST data website. You will need to be connected to the Internet.
# DataDeps is a utility for downloading data on demand from a URL
using DataDeps
# Where do I want to put the data?
ENV["DATADEPS_LOAD_PATH"] = joinpath(datadir(), "exp_raw")
ENV["DATADEPS_ALWAYS_ACCEPT"]="true"
# Register the data sets using the names "MnDeepSeaNodule" and "MnDeepSeaNoduleStandards"
register(DataDep("MnDeepSeaNodule",
"""
Dataset: Deep Sea Manganese Nodule
Acquired by: Nicholas W.M. Ritchie
License: CC-SA 3.0
""",
raw"https://data.nist.gov/od/ds/mds2-2467/MnNodule.tar.gz",
"5b5b6623b8f4daca3ff3073708442ac5702ff690aa12668659875ec5642b458d",
post_fetch_method=DataDeps.unpack
));
register(DataDep("MnDeepSeaNoduleStandards",
"""
Dataset: Standard Spectra for Deep Sea Manganese Nodule Dataset
Acquired by: Nicholas W.M. Ritchie
License: CC-SA 3.0
""",
raw"https://data.nist.gov/od/ds/mds2-2467/MnNodule_Standards.tar.gz",
"69283ba72146932ba451e679cf02fbd6b350f96f6d012d50f589ed9dd2e35f1a",
post_fetch_method=DataDeps.unpack
));
The data is in RPL/RAW format - a format designed by David Bright at NIST and available from many EDS vendor's software. In addition to the data in the RPL/RAW format, we also need to provide spectrum meta-data and the detector energy scale.
raw = readrplraw(joinpath(datadep"MnDeepSeaNodule","map[15]"))
les = LinearEnergyScale(0.0, 10.0)
props = Dict{Symbol,Any}(
:LiveTime => 0.72*4.0*18.0*3600.0/(1024.0*1024.0),
:BeamEnergy => 20.0e3,
:TakeOffAngle => deg2rad(35.0),
:ProbeCurrent => 1.0,
:Name => "Deep Sea Mn Nodule",
)
hs = HyperSpectrum(les, props, raw, fov = (0.2, 0.2))
1024 × 1024 HyperSpectrum{UInt16,(:Y, :X)}[Deep Sea Mn Nodule), 0.0 + 10.0⋅ch eV, 2048 ch]
It is easy to generate crude elemental maps from HyperSpectrum structs. These are simple ROI maps constructed by summing a range of channels around the central channel for the specified X-ray transition.
hs[n"Mn K-L3"]
Much like spectra, hyperspectra have properties. However, the same set of properties apply to all spectra within a hyperspectrum. There is one exception, :LiveTime, which can be an array of the same dimensionality as the hyperspectrum to represent a unique live-time per pixel.
NeXLCore.properties(hs)
Dict{Symbol, Any} with 5 entries:
:Name => "Deep Sea Mn Nodule"
:BeamEnergy => 20000.0
:LiveTime => 0.177979
:TakeOffAngle => 0.610865
:ProbeCurrent => 1.0
# This speeds up the processing by binning a 1024 x 1024 hyperspectrum to 1024/block_size x 1024/block_size.
block_size = 1
hs= block_size > 1 ? block(hs, (block_size, block_size)) : hs
1024 × 1024 HyperSpectrum{UInt16,(:Y, :X)}[Deep Sea Mn Nodule), 0.0 + 10.0⋅ch eV, 2048 ch]
Unfortunately, ROI maps suffer from many shortcomings including
While they are convenient as a first glimse at the data, they can be very deceptive and should rarely be included in reports.
Instead, it is better to fit the point spectra within the hyperspectrum using measured references.
hs_elms = [ n"Ag", n"Al", n"Ba", n"C", n"Ca", n"Ce", n"Cr", n"Cu", n"Fe", n"S", n"P", n"K", n"Mg", n"O", n"Mn", n"Na", n"Cl", n"Ni", n"Si", n"Ti", n"Zn" ]
plot(maxpixel(hs), klms = hs_elms, xmax=10.0e3)
print(symbol.(sort(hs_elms)))
["C", "O", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "K", "Ca", "Ti", "Cr", "Mn", "Fe", "Ni", "Cu", "Zn", "Ag", "Ba", "Ce"]
refs = references( [
reference(n"C", joinpath(datadep"MnDeepSeaNoduleStandards", "C std.msa")),
reference(n"O", joinpath(datadep"MnDeepSeaNoduleStandards", "MgO std.msa")),
reference(n"Na", joinpath(datadep"MnDeepSeaNoduleStandards", "NaCl std.msa")),
reference(n"Mg", joinpath(datadep"MnDeepSeaNoduleStandards", "Mg std.msa")),
reference(n"Al", joinpath(datadep"MnDeepSeaNoduleStandards", "Al std.msa")),
reference(n"Si", joinpath(datadep"MnDeepSeaNoduleStandards", "Si std.msa")),
reference(n"P", joinpath(datadep"MnDeepSeaNoduleStandards", "GaP std.msa")),
reference(n"S", joinpath(datadep"MnDeepSeaNoduleStandards", "FeS2 std.msa")),
reference(n"Cl", joinpath(datadep"MnDeepSeaNoduleStandards", "NaCl std.msa")),
reference(n"K", joinpath(datadep"MnDeepSeaNoduleStandards", "KBr std.msa")),
reference(n"Ca", joinpath(datadep"MnDeepSeaNoduleStandards", "CaF2 std.msa")),
reference(n"Ti", joinpath(datadep"MnDeepSeaNoduleStandards", "Ti std.msa")),
reference(n"Cr", joinpath(datadep"MnDeepSeaNoduleStandards", "Cr std.msa")),
reference(n"Mn", joinpath(datadep"MnDeepSeaNoduleStandards", "Mn std.msa")),
reference(n"Fe", joinpath(datadep"MnDeepSeaNoduleStandards", "Fe std.msa")),
reference(n"Ni", joinpath(datadep"MnDeepSeaNoduleStandards", "Ni std.msa")),
reference(n"Cu", joinpath(datadep"MnDeepSeaNoduleStandards", "Cu std.msa")),
reference(n"Zn", joinpath(datadep"MnDeepSeaNoduleStandards", "Zn std.msa")),
reference(n"Ag", joinpath(datadep"MnDeepSeaNoduleStandards", "Ag std.msa")),
reference(n"Ba", joinpath(datadep"MnDeepSeaNoduleStandards", "BaF2 std.msa"))
], 132.0)
asa(DataFrame, refs)
[ Info: A material containing Ba and F is not a suitable reference for "Ba M5-N3 + 16 others" due to a peak interference.
| Row | Spectrum | Beam Energy (keV) | Probe Current (nA) | Live Time (s) | Material | Lines | ROI | P-to-B | S-to-N |
|---|---|---|---|---|---|---|---|---|---|
| SubStrin… | Float64 | Float64 | Float64 | Material… | Array… | UnitRang… | Float64 | Float64 | |
| 1 | C std | 20.0 | 0.97372 | 1184.6 | C[C=1.0000,1.90 g/cm³] | C K-L3 + 1 other | 16:40 | 24.143 | 15279.7 |
| 2 | MgO std | 20.0 | 0.9748 | 1161.48 | MgO[Mg=0.6030,O=0.3970,5.00 g/cm³] | O K-L3 + 1 other | 39:66 | 12.7889 | 8706.54 |
| 3 | NaCl std | 20.0 | 0.97454 | 1156.51 | NaCl[Na=0.3934,Cl=0.6066,5.00 g/cm³] | Na K-L3 + 1 other | 89:120 | 6.21839 | 9555.21 |
| 4 | Mg std | 20.0 | 0.97475 | 1149.05 | Mg[Mg=1.0000,1.74 g/cm³] | Mg K-L3 + 1 other | 110:142 | 17.4471 | 34147.0 |
| 5 | Al std | 20.0 | 0.97396 | 1151.25 | Al[Al=1.0000,1.00 g/cm³] | Al K-L3 + 3 others | 132:169 | 19.6389 | 36840.4 |
| 6 | Si std | 20.0 | 0.97465 | 1149.58 | Si[Si=1.0000,2.95 g/cm³] | Si K-L3 + 3 others | 157:198 | 25.0393 | 40461.5 |
| 7 | GaP std | 20.0 | 0.97451 | 1161.35 | GaP[P=0.3076,Ga=0.6924,5.00 g/cm³] | P K-L3 + 3 others | 184:230 | 8.27172 | 9393.98 |
| 8 | FeS2 std | 20.0 | 0.97334 | 1160.73 | FeS2[S=0.5345,Fe=0.4655,5.00 g/cm³] | S K-L3 + 3 others | 212:264 | 16.8388 | 21532.6 |
| 9 | NaCl std | 20.0 | 0.97454 | 1156.51 | NaCl[Na=0.3934,Cl=0.6066,5.00 g/cm³] | Cl K-L3 + 3 others | 243:300 | 23.512 | 25791.1 |
| 10 | KBr std | 20.0 | 0.97392 | 1147.26 | KBr[Br=0.6714,K=0.3286,5.00 g/cm³] | K K-L3 + 3 others | 310:379 | 7.04575 | 9501.98 |
| 11 | CaF2 std | 20.0 | 0.97401 | 1167.54 | CaF2[Ca=0.5133,F=0.4867,3.18 g/cm³] | Ca K-L3 + 3 others | 347:422 | 17.5274 | 19125.9 |
| 12 | Ti std | 20.0 | 0.9752 | 1160.26 | Ti[Ti=1.0000,5.00 g/cm³] | Ti L3-M5 + 13 others | 28:65 | 4.35094 | 3026.39 |
| 13 | Ti std | 20.0 | 0.9752 | 1160.26 | Ti[Ti=1.0000,5.00 g/cm³] | Ti K-L3 + 3 others | 427:516 | 20.5636 | 26561.6 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 22 | Fe std | 20.0 | 0.97436 | 1164.66 | Fe[Fe=1.0000,7.87 g/cm³] | Fe K-M3 + 3 others | 680:732 | 3.77854 | 3010.53 |
| 23 | Ni std | 20.0 | 0.97258 | 1164.08 | Ni[Ni=1.0000,5.00 g/cm³] | Ni L3-M5 + 13 others | 62:108 | 10.0052 | 11386.8 |
| 24 | Ni std | 20.0 | 0.97258 | 1164.08 | Ni[Ni=1.0000,5.00 g/cm³] | Ni K-L3 + 1 other | 718:778 | 16.0413 | 14413.1 |
| 25 | Ni std | 20.0 | 0.97258 | 1164.08 | Ni[Ni=1.0000,5.00 g/cm³] | Ni K-M3 + 3 others | 799:855 | 3.18328 | 2401.56 |
| 26 | Cu std | 20.0 | 0.97242 | 1159.85 | Cu[Cu=1.0000,8.90 g/cm³] | Cu L3-M5 + 13 others | 69:117 | 15.3703 | 19251.1 |
| 27 | Cu std | 20.0 | 0.97242 | 1159.85 | Cu[Cu=1.0000,8.90 g/cm³] | Cu K-L3 + 1 other | 773:836 | 14.0936 | 12296.1 |
| 28 | Cu std | 20.0 | 0.97242 | 1159.85 | Cu[Cu=1.0000,8.90 g/cm³] | Cu K-M3 + 3 others | 862:920 | 2.93067 | 2123.89 |
| 29 | Zn std | 20.0 | 0.97292 | 1157.03 | Zn[Zn=1.0000,5.00 g/cm³] | Zn L3-M5 + 13 others | 76:126 | 17.3497 | 23748.8 |
| 30 | Zn std | 20.0 | 0.97292 | 1157.03 | Zn[Zn=1.0000,5.00 g/cm³] | Zn K-L3 + 1 other | 831:896 | 12.5069 | 10519.0 |
| 31 | Zn std | 20.0 | 0.97292 | 1157.03 | Zn[Zn=1.0000,5.00 g/cm³] | Zn K-M3 + 3 others | 928:988 | 2.67558 | 1862.94 |
| 32 | Ag std | 20.0 | 0.98647 | 921.028 | Ag[Ag=1.0000,5.00 g/cm³] | Ag L3-M5 + 25 others | 247:393 | 5.28945 | 11615.2 |
| 33 | BaF2 std | 20.0 | 0.98986 | 1171.3 | BaF2[Ba=0.7833,F=0.2167,5.00 g/cm³] | Ba L3-M5 + 29 others | 377:617 | 2.85456 | 7169.28 |
Now we will fit the references to the hyperspectrum unknown. This is simple enough:
kratios_fast = fit_spectrum(hs, refs, mode = :Fast)
This is time consuming so once we've done it, we want to store the result in a file and reload the result rather than recalculate it. DrWatsons provides a mechanism to do this based on a function produce_or_load(....). You specify where to put the data and it builds a filename based on the prefix, suffix and the configuration paramters in @dict(mode, x_dim, y_dim). The suffix specifies that the data is written using the JLD2 library in HDF5 format.
# kratios_fast = @time fit_spectrum(hs, refs, mode = :Fast) # Single threaded took 92 seconds
using FileIO, Parameters
mode = :Fast
x_dim, y_dim = size(hs)
data, file = produce_or_load(
joinpath(datadir(),"exp_pro"),
@dict(mode, x_dim, y_dim),
prefix="kratios", suffix="jld2"
) do config
Dict("kratios" => @time fit_spectrum(hs, refs, mode = mode))
end
@show file
kratios_fast = data["kratios"]
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\kratios_mode=Fast_x_dim=1024_y_dim=1024.jld2"
33-element Vector{KRatios}:
k[C, C K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[MgO, O K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[NaCl, Na K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[Mg, Mg K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[Al, Al K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[Si, Si K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[GaP, P K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[FeS2, S K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[NaCl, Cl K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[KBr, K K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[CaF2, Ca K-L3 + 3 others] = Float64[ (1024, 1024) ]
k[Ti, Ti L3-M5 + 13 others] = Float64[ (1024, 1024) ]
k[Ti, Ti K-L3 + 3 others] = Float64[ (1024, 1024) ]
⋮
k[Fe, Fe K-M3 + 3 others] = Float64[ (1024, 1024) ]
k[Ni, Ni L3-M5 + 13 others] = Float64[ (1024, 1024) ]
k[Ni, Ni K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[Ni, Ni K-M3 + 3 others] = Float64[ (1024, 1024) ]
k[Cu, Cu L3-M5 + 13 others] = Float64[ (1024, 1024) ]
k[Cu, Cu K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[Cu, Cu K-M3 + 3 others] = Float64[ (1024, 1024) ]
k[Zn, Zn L3-M5 + 13 others] = Float64[ (1024, 1024) ]
k[Zn, Zn K-L3 + 1 other] = Float64[ (1024, 1024) ]
k[Zn, Zn K-M3 + 3 others] = Float64[ (1024, 1024) ]
k[Ag, Ag L3-M5 + 25 others] = Float64[ (1024, 1024) ]
k[BaF2, Ba L3-M5 + 29 others] = Float64[ (1024, 1024) ]
The result of fit_spectrum(...) is a set k-ratio maps. Which can be readily visualized and represent a more quantitative perspective on the spectrum data. The principle advantage being that the data is background corrected so that continuum is not mistaken for a small quantity of an element.
To visualize the k-ratios, it is best to pick one k-ratio per element using optimizeks(...) and then normalize the k-ratios on a point-by-point basis using normalizek(...). This ensures that all the k-ratios are between 0 and 1. The function asa(ElementalMap, ...) facilitates this process by converting a Vector{KRatios} into a dictionary of images indexed by Element instances.
using Colors
imgs = asa(ElementalMap, kratios_fast, Gray)
imgs[n"Mn"]
We can perform matrix correction on the Vector{KRatios} using the quantify(...) method. It returns a Materials struct which is a memory efficient way of representing a multi-dimensional array of Material.
# quant_fast=quantify(kratios_fast, name=hs[:Name], ty=Float32) # Takes ~80 minute for 1024 x 1024 (~5 ms per pixel)
data, file = produce_or_load(
joinpath(datadir(),"exp_pro"),
@dict(mode, x_dim, y_dim),
prefix="quantify", suffix="jld2"
) do config
Dict("mass_fractions" => @time quantify(kratios_fast, name=hs[:Name], ty=Float32))
end
@show file
quant_fast = data["mass_fractions"]
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\quantify_mode=Fast_x_dim=1024_y_dim=1024.jld2"
The Materials struct can be indexed at a specific coordinate to extract the composition at that coordinate as a Material.
quant_fast[128,384]
Deep Sea Mn Nodule[128,384][Ba=0.0008,Ca=0.0087,Mg=0.0006,O=0.2713,Cl=0.0064,K=0.0008,Ti=0.0011,Ag=0.0017,C=1.0478,Na=0.0054,Si=0.0414,Cu=0.0088,Mn=0.0039,Al=0.0143,Fe=0.0092]
Using broadcast, we can apply Material functions to Materials objects. The results is an Array with the same shape but with the result contents.
at_ms=analyticaltotal.(quant_fast)
1024×1024 Matrix{Float64}:
1.12114 0.931418 1.0184 1.12544 … 0.798097 0.798483 0.750225
0.877755 0.939474 0.87016 0.839677 0.826642 0.753609 0.814797
0.944509 0.982147 0.936352 0.894357 0.789687 0.802227 0.806136
0.877553 0.834916 0.877479 0.952697 0.69093 0.722851 0.736526
0.83982 0.831247 0.790591 0.944809 0.604645 0.658924 0.682592
0.874638 0.908482 0.92161 0.851666 … 0.673269 0.714052 0.675014
0.920448 0.973551 0.997805 0.919651 0.822391 0.687747 0.698407
0.901042 0.892611 0.878813 0.943202 0.794712 0.749333 0.647348
0.840059 0.901996 0.919389 0.944544 0.783919 0.788288 0.667246
0.904216 0.947249 0.971246 0.89693 0.926894 0.825388 0.742284
0.972421 1.02393 0.94819 0.948975 … 0.862409 0.826175 0.766345
0.932931 0.945948 0.882228 0.975659 0.782198 0.80973 0.804723
0.916943 0.905414 0.837847 0.824939 0.891919 0.883625 0.82066
⋮ ⋱
0.882333 0.928941 0.937709 0.994533 0.965608 0.838502 0.799263
1.00358 0.993153 0.931432 0.955147 0.736347 0.883015 0.900127
1.02583 1.06478 0.964631 1.00421 0.840045 0.864325 0.949038
1.03194 0.947343 0.85004 1.06743 … 0.912345 0.97687 0.922914
1.03552 1.02994 1.01893 1.00077 0.987512 1.06054 1.09555
1.0135 1.03027 0.980763 0.945134 0.809482 0.875999 1.08297
1.10909 1.11128 1.03744 0.993881 0.602245 0.734092 0.860694
1.15768 1.01753 1.01024 1.04317 0.517874 0.572694 0.701184
0.989872 1.1088 1.00917 1.01493 … 0.520139 0.556505 0.810547
1.0812 1.00858 0.999245 1.00088 0.606844 0.731831 0.94359
0.911391 0.961544 0.924331 0.962723 0.745007 0.92867 0.979702
0.982772 1.04914 1.00912 0.988079 0.78446 0.890835 0.926269
using Statistics
mean(at_ms), std(at_ms)
(1.0447547161174755, 0.1606245028098845)
To convert the mass-fractions to images, we need to ensure all mass-fractions are on the range [0,1]. We can do this by normalizing the individual points to an analytical total of unity.
quant_fast_norm=asnormalized(quant_fast)
Now let's visualize this. In Julia, it is trivial to convert matrices representing values on the range [0,1] to images.
We will extract the plane of data corresponding to an element by indexing it with an Element. Log3Band is a function that converts numbers between 0 and 1 to RGB values. The mapping is displayed below in the legend.
Log3Band.(quant_fast_norm[n"Ni"])
loadlegend("Log3BandBright.png")
The Ni elemental map can be seen to represent primarily values ranging from a few percent down to zero. Albeit, this spectrum image took 18 hours to collect but the depth of information might surprise many who are used to thinking of spectrum images as a crude tool.
Let's output various perspectives from each element to the 'plots' directory. The Gray function converts values on the range [0,1] to gray-scale values.
# Fast k-ratio maps
let imgs = asa(ElementalMap, kratios_fast, Log3Band)
for (el, img) in imgs
save(joinpath(plotsdir(),"$(symbol(el)) - Fast Log3Band.png"), img)
end
end
# Full k-ratio maps
let imgs = asa(ElementalMap, kratios_fast, Gray)
for (el, img) in imgs
save(joinpath(plotsdir(),"$(symbol(el)) - Fast Gray.png"), img)
end
end
# Full quantitative maps - Fast fit
begin
for el in keys(quant_fast_norm)
save(joinpath(plotsdir(),"$(symbol(el)) - Quant Fast Gray.png"), Gray.(quant_fast_norm[el]))
save(joinpath(plotsdir(),"$(symbol(el)) - Quant Fast Log3Band.png"), Log3Band.(quant_fast_norm[el]))
end
end
All of these results have been based on the mode=:Fast fit of the standard spectra to the hyperspectrum data.
Now we will perform a mode=:Full fit of the HyperSpectrum data. This takes a couple of orders of magnitude more time than the mode=:Fast fit. However, the :Full mode won't return negative k-ratios and, in general, performs a much more careful fit.
# kratios_full = @time fit_spectrum(hs, refs, mode = :Full) # Single threaded took 6132 seconds
mode = :Full
data, file = produce_or_load(
joinpath(datadir(),"exp_pro"),
@dict(mode, x_dim, y_dim),
prefix="kratios", suffix="jld2"
) do config
Dict("kratios" => @time fit_spectrum(hs, refs, mode = mode))
end
@show file
kratios_full = data["kratios"]
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\kratios_mode=Full_x_dim=1024_y_dim=1024.jld2"
33-element Vector{KRatios}:
k[C, C K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[MgO, O K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[NaCl, Na K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[Mg, Mg K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[Al, Al K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[Si, Si K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[GaP, P K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[FeS2, S K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[NaCl, Cl K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[KBr, K K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[CaF2, Ca K-L3 + 3 others] = Float32[ (1024, 1024) ]
k[Ti, Ti L3-M5 + 13 others] = Float32[ (1024, 1024) ]
k[Ti, Ti K-L3 + 3 others] = Float32[ (1024, 1024) ]
⋮
k[Fe, Fe K-M3 + 3 others] = Float32[ (1024, 1024) ]
k[Ni, Ni L3-M5 + 13 others] = Float32[ (1024, 1024) ]
k[Ni, Ni K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[Ni, Ni K-M3 + 3 others] = Float32[ (1024, 1024) ]
k[Cu, Cu L3-M5 + 13 others] = Float32[ (1024, 1024) ]
k[Cu, Cu K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[Cu, Cu K-M3 + 3 others] = Float32[ (1024, 1024) ]
k[Zn, Zn L3-M5 + 13 others] = Float32[ (1024, 1024) ]
k[Zn, Zn K-L3 + 1 other] = Float32[ (1024, 1024) ]
k[Zn, Zn K-M3 + 3 others] = Float32[ (1024, 1024) ]
k[Ag, Ag L3-M5 + 25 others] = Float32[ (1024, 1024) ]
k[BaF2, Ba L3-M5 + 29 others] = Float32[ (1024, 1024) ]
Let's add the :Full fit k-ratio data images to the 'plots' directory for comparison.
using Colors
let imgs = asa(ElementalMap, kratios_full, Log3Band)
for (el, img) in imgs
save(joinpath(plotsdir(),"$(symbol(el)) - Full Log3Band.png"), img)
end
end
let imgs = asa(ElementalMap, kratios_full, Gray)
for (el, img) in imgs
save(joinpath(plotsdir(),"$(symbol(el)) - Full Gray.png"), img)
end
end
Let's compare the results from the full fit with the results from the fast fit.
imgs_full = asa(ElementalMap, kratios_full, Log3Band)
els = sort([ keys(imgs_full)...])
print(symbol.(els))
["C", "O", "Na", "Mg", "Al", "Si", "P", "S", "Cl", "K", "Ca", "Ti", "Cr", "Mn", "Fe", "Ni", "Cu", "Zn", "Ag", "Ba"]
map( el->imgs_full[el], els)
begin
imgs_fast = asa(ElementalMap, kratios_fast, Log3Band)
map( el->imgs_fast[el], els)
end
Now let's quantify the :Full fit data. Performing the matrix correction is actually much slower than performing the :Full fit.
Again, the results are a Materials struct.
# quant_full = @time quantify(kratios_full, name=hs[:Name])
data, file = produce_or_load(
joinpath(datadir(),"exp_pro"),
@dict(mode, x_dim, y_dim),
prefix="quantify", suffix="jld2"
) do config
Dict("mass_fractions" => @time quantify(kratios_full, name=hs[:Name], ty=Float32))
end
@show file
quant_full = data["mass_fractions"]
[ Info: File c:\Users\nritchie\Desktop\IUMAS\data\exp_pro\quantify_mode=Full_x_dim=1024_y_dim=1024.jld2 does not exist. Producing it now...
4467.479611 seconds (32.22 G allocations: 7.521 TiB, 16.99% gc time, 0.03% compilation time: 49% of which was recompilation)
┌ Warning: The Git repository ('c:\Users\nritchie\Desktop\IUMAS') is dirty! Appending -dirty to the commit ID. └ @ DrWatson C:\Users\nritchie\.julia\packages\DrWatson\JACNB\src\saving_tools.jl:64
file = "c:\\Users\\nritchie\\Desktop\\IUMAS\\data\\exp_pro\\quantify_mode=Full_x_dim=1024_y_dim=1024.jld2"
[ Info: File c:\Users\nritchie\Desktop\IUMAS\data\exp_pro\quantify_mode=Full_x_dim=1024_y_dim=1024.jld2 saved.
map(el->Log3Band.(quant_full[el]), els)
Clearly, the light element maps are where we see the influence of matrix correction. When X-ray absorption is small, the k-ratio is generally a good approximation of the mass-fraction. However, for low energy X-rays, the absorption can be strong leading to a significant underestimation of the mass fraction. We see this particularly in the carbon data.
Now write these images to disk.
begin
quant_full_norm=asnormalized(quant_full)
for el in keys(quant_full_norm)
save(joinpath(plotsdir(),"$(symbol(el)) - Quant Full Gray.png"), Gray.(quant_full_norm[el]))
save(joinpath(plotsdir(),"$(symbol(el)) - Quant Full Log3Band.png"), Log3Band.(quant_full_norm[el]))
end
end
quant_full[n"Mn"]
1024×1024 Matrix{Float32}:
0.223394 0.311174 0.332411 … 0.192086 0.18319 0.182768
0.269441 0.277169 0.307894 0.187049 0.187401 0.198357
0.27579 0.272186 0.242913 0.1994 0.215164 0.195631
0.284984 0.254272 0.251819 0.221909 0.19468 0.228718
0.271912 0.248145 0.231546 0.181427 0.315283 0.155952
0.00486168 0.235095 0.232646 … 0.209557 0.0161902 0.198351
0.216475 0.192771 0.306901 0.200423 0.165684 0.218023
0.220113 0.231561 0.285585 0.213188 0.22912 0.21971
0.248034 0.252878 0.272362 0.19302 0.214149 0.281862
0.24977 0.260708 0.239762 0.301234 0.185244 0.209922
0.266883 0.271193 0.297616 … 0.165329 0.190385 0.211984
0.301004 0.146325 0.274044 0.26669 0.218611 0.207211
0.29647 0.303459 0.285447 0.265315 0.248976 0.217168
⋮ ⋱
0.324508 0.324034 0.26046 0.00610759 0.00819175 0.0323414
0.223967 0.298605 0.231759 0.247855 0.0136075 0.0366124
0.267053 0.2395 0.00541237 0.00076621 0.0255185 0.0395083
0.0377666 0.161325 0.092179 … 0.0 0.416347 0.00254732
0.147561 0.097022 0.0636286 0.193515 0.00693751 0.00533547
0.0941413 0.0659752 0.0677662 0.00869319 0.00985284 0.00891134
0.101526 0.07201 0.104355 0.0264528 0.017331 0.00477545
0.27833 0.129675 0.188148 0.036561 0.0274113 0.0151129
0.15678 0.207637 0.274506 … 0.0463395 0.0424375 0.334446
0.225777 0.146953 0.385326 0.069456 0.0334889 0.0160826
0.208606 0.294549 0.346879 0.0603972 0.0310457 0.0142638
0.220558 0.30014 0.321251 0.0684462 0.0361965 0.0130948
using Statistics
mean(analyticaltotal.(quant_full)), mean(analyticaltotal.(quant_fast))
(1.0397453436659276, 1.0447547161174755)
To explore the data, let's plot histograms of the mass fractions from each element.
colors = distinguishable_colors(length(els), colorant"light blue", lchoices=range(50, stop=100, length=15))
plot(
( layer(x=quant_full[el], Geom.histogram, Theme(default_color=colors[i]), alpha=[0.6]) for (i, el) in enumerate(els) )...,
Guide.xlabel("Mass Fraction"), Guide.ylabel("Counts"),
Guide.manual_color_key("Element", symbol.(els), colors),
Coord.cartesian(xmin=0.0, xmax=1.0, ymax=16.0e4)
)
plot(
( layer(x=quant_full[el], Geom.histogram, Theme(default_color=colors[i]), alpha=[0.6]) for (i, el) in enumerate(els) )...,
Guide.xlabel("Mass Fraction"), Guide.ylabel("Counts"),
Guide.manual_color_key("Element", symbol.(els), colors),
Coord.cartesian(xmin=0.0, xmax=0.05, ymax=16.0e4)
)
We can easily query the quantified data using broadcast syntax to create BitMatrix based on a boolean comparison.
To determine which pixels have a mass-fraction of Si > 5% we create a mask. In this BitMatrix mask, 1 represents a true comparison and 0 a false comparison.
mask_si = quant_full[n"Si"] .> 0.05
1024×1024 BitMatrix: 0 0 0 1 0 1 1 1 1 0 0 0 0 … 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 … 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 1 … 0 0 0 0 0 0 0 1 1 1 0 0 0 1 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 0 0 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 1 1 0 0 1 1 1 1 1 1 1 1 1 0 0 1 0 1 0 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 1 0 1 1 0 … 1 0 0 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 1 0 1 1 1 1 1 0 1 1 1 1 1 0 0 1 1 1 1 1 0 0 0 1 1 1 1 0 1 1 0 1 1 1 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 … 1 1 1 1 1 1 1 1 0 0 0 1 0 1 1 0 0 0 0 0 0 0 0 0 0 0 1 0 1 1 0 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 1 1 1 1 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 1 1 1
Using similar syntax we can determine how many pixels have at least 50%, 10%, 5% and 1% of an element?
DataFrame(
Element=symbol.(els),
P50=map(el->count(quant_full[el] .> 0.5),els),
P10=map(el->count(quant_full[el] .> 0.1),els),
P5=map(el->count(quant_full[el] .> 0.05),els),
P1=map(el->count(quant_full[el] .> 0.01),els)
)
| Row | Element | P50 | P10 | P5 | P1 |
|---|---|---|---|---|---|
| String | Int64 | Int64 | Int64 | Int64 | |
| 1 | C | 86704 | 720858 | 1019909 | 1045263 |
| 2 | O | 35574 | 1034995 | 1040636 | 1045908 |
| 3 | Na | 0 | 39 | 5515 | 603161 |
| 4 | Mg | 0 | 96 | 896 | 727196 |
| 5 | Al | 0 | 5380 | 39062 | 742418 |
| 6 | Si | 5 | 121297 | 270774 | 942104 |
| 7 | P | 0 | 213 | 653 | 2336 |
| 8 | S | 0 | 1 | 28 | 115253 |
| 9 | Cl | 0 | 0 | 1 | 23493 |
| 10 | K | 0 | 453 | 3768 | 269937 |
| 11 | Ca | 0 | 1190 | 5723 | 809920 |
| 12 | Ti | 0 | 159 | 469 | 4871 |
| 13 | Cr | 0 | 0 | 3 | 121 |
| 14 | Mn | 184 | 863934 | 918021 | 990386 |
| 15 | Fe | 33 | 182310 | 506314 | 960971 |
| 16 | Ni | 2 | 10 | 15 | 543173 |
| 17 | Cu | 0 | 6 | 11 | 416359 |
| 18 | Zn | 2 | 75 | 214 | 85702 |
| 19 | Ag | 0 | 1 | 12 | 17231 |
| 20 | Ba | 0 | 212 | 1064 | 67902 |
A mask can then be applied to extract and sum the requisite spectra within the hyperspectrum.
plot( sum(hs, quant_full[n"Si"] .> 0.05), klms=els, xmax=10.0e3)
Let's do something similar for Cu > 10%.
plot( sum(hs, quant_full[n"Cu"] .> 0.1), klms=els, xmax=10.0e3)
This is just a flavor of what you can do with spectra and hyperspectra using Julia, the Julia library infrastructure including the NeXL libraries. There are many machine learning, chemometric, image analysis and other libraries that can be readily applied.